/* Include files */
#include "math.h"
#include "mex.h"


/* Checks whether inputs are matrices */
#define IS_REAL_2D_FULL_DOUBLE(P) (!mxIsComplex(P) && \
mxGetNumberOfDimensions(P) == 2 && !mxIsSparse(P) && mxIsDouble(P))
#define IS_REAL_SCALAR(P) (IS_REAL_2D_FULL_DOUBLE(P) && mxGetNumberOfElements(P) == 1)

/*************************************************/
/* MAIN                                          */
/*************************************************/

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
	/* Define output matrix */
	#define Share    	plhs[0]

	/* Define input matrices */
	#define Firm 		prhs[0]
	#define Market      prhs[1]
	#define J 			prhs[2]

	int i, j, m;

	/* Check the number of inputs is correct (Aoffer,Baccept,QbarA) */
	if(nrhs < 3)
		mexErrMsgTxt("Too few input arguments.");
	else if(nrhs > 3)
		mexErrMsgTxt("Too many input arguments.");

	if(!IS_REAL_2D_FULL_DOUBLE(Firm))
		mexErrMsgTxt("Firm must be real 2D full double array.");
	if(!IS_REAL_2D_FULL_DOUBLE(Market))
		mexErrMsgTxt("Market must be real 2D full double array.");
	if(!IS_REAL_2D_FULL_DOUBLE(J))
		mexErrMsgTxt("J must be real 2D full double array.");

	/* Dimensions and pointers for inputs */
	int N = mxGetM(Firm);
	int M = mxGetM(Market);

	double *Firm_pt 	= mxGetPr(Firm);
	double *Market_pt 	= mxGetPr(Market);
	double *J_pt 		= mxGetPr(J);

	int Jm = J_pt[0]*M;
	int Jn = J_pt[0]*N;

	/* Output matrix */
	Share = mxCreateDoubleMatrix(Jm,1,mxREAL);
	double *Share_pt = mxGetPr(Share);
	for (j = 0; j < Jm; j++)
	{
		Share_pt[j] = 0;
	}

/***************************************************************/
/* Begin Filling       */
/***************************************************************/

	int shareID = 0;
	int beg = 0;
	int end = 0;
	int denom = 0;

	for(m = 0; m < M; m++)
	{
		beg = Market_pt[m + M*1];
		end = Market_pt[m + M*2];
		for (j = beg - 1; j < end; j++)
		{			
			if (Firm_pt[j] > 0)
			{
				shareID = Firm_pt[j] + m*J_pt[0] - 1;
				denom = Market_pt[m + M*0];
				Share_pt[shareID] = Share_pt[shareID] + 1/(float)denom;
			}
		}
	}

}
